Unilateral interactions in granular packings: A model for the anisotropy modulus 
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Unilateral interparticle interactions have an effect on the elastic response of granular materials 
due to the opening and closing of contacts during quasi-static shear deformations. A simplified 
model is presented, for which constitutive relations can be derived. For biaxial deformations the 
elastic behavior in this model involves three independent elastic moduli: bulk, shear, and anisotropy 
modulus. The bulk and the shear modulus, when scaled by the contact density, are independent of 
the deformation. However, the magnitude of the anisotropy modulus is proportional to the ratio 
between shear and volumetric strain. Sufficiently far from the jamming transition, when corrections 
due to non-affine motion become weak, the theoretical predictions are qualitatively in agreement 
with simulation results. 

PACS numbers: 45.70.-n, 46.25.-y, 83.80.Fg 
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I. INTRODUCTION 

Understanding the mechanical response of granular 
materials is still one of the remaining challenges in ma- 
terials science and physics. The research is motivated 
by many industrial and geophysical applications [H, 
A granular packing at rest does not behave like an ordi- 
nary elastic solid, because the relation between stress and 
strain is nonlinear, depends on the fabric and on the load- 
ing path, and gives rise to energy dissipation [sUEIl . These 
properties have been investigated for packings of spheres, 
where they could be traced back to the nonlinearity of 
Hertzian contacts, disorder, and Coulomb friction law. 
In these studies the contacts between the spheres were 
supposed to remain closed under the strain. 

However, there is another source of the nonlinear elas- 
tic response in granular media, which has not been stud- 
ied in such detail yet @. Due to the absence of attrac- 
tive contact forces in dry granular media, the contacts 
can open, and thus do not transmit any elastic restoring 
nor frictional force. Even in the case of linear (Hookean) 
elasticity on the particle level (e.g. for packings of paral- 
lel cylinders), the vanishing of the elastic response under 
tension renders the macroscopic elastic behavior nonlin- 
ear. It is this nonlinearity which we analyze in this paper. 

It is known that the contact and force networks in 
sheared granular materials are anisotropic 043- The 
present analytical approach links the fabric anisotropy 
to the shear deformation for fixed volumetric strain. 

In the following, we investigate the incremental lin- 
ear response of a pre-strained two dimensional packing 
of disks and relate the elements of the stiffness tensor to 
the mean packing properties and the probability distribu- 
tion of contact orientations. Then, analytical expressions 
for the elastic moduli of isotropic and sheared contact 
networks are derived, and the results are compared with 
numerical simulations. 



II. LINEAR ELASTIC RESPONSE 

The elastic response of solids to deformations is classi- 
cally expressed in terms of the relationship between the 
stress tensor <t and the strain tensor e. Let us assume 
that the relation between the incremental stress tensor 
Scr and the incremental strain tensor Se can be expressed 
by analytical functions F. . as 
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where F^^ satisfy (Je,, =0, V fc, / e {1, d})=0. Here, 
d is the dimension of the system. In the limit of small 
deformations, Eq. ([T|) can be approximated by Taylor 
expansion to first order in the strain increment 
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where C is a tensor of rank 4, referred to as the stiff- 
ness or elasticity tensor. Equation ([2]) can be regarded 
as the generalization of Hooke's law that describes the 
linear response of the medium to external perturbations. 
Although C has 16 elements in a two-dimensional system, 
symmetry considerations for the stress and strain tensors 
imply that there are usually less independent elements. 
For example, the symmetry of a completely isotropic ma- 
terial requires that C has only two independent elements, 
commonly represented by the Lame coefficients or, alter- 
natively, by Poisson's ratio and Young's modulus. 

In this section, it is recalled, how an approximate ex- 
pression for the stiffness tensor of a dense assembly of 
grains can be obtained, which allows to calculate the elas- 
tic moduli from mean packing properties and the proba- 
bility distribution of contact orientations [l^, [HI ■ 

The schematic figure [IJa) shows the undeformed 
shapes of two particles at positions A and B. Their over- 
lap is a measure of the elastic deformation at their con- 
tact c. The normal and tangential contact unit vectors 
are denoted by n"^ and f^, and the branch vector l'^ con- 
nects the centers of the particles. Each contact force is 
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FIG. 1: (a) Geometry of the contact c between the particles A 
and B including the branch vector l'^, the normal unit vector 
n'^, and the tangential unit vector f^. (b) n'^ and in an 
arbitrarily chosen Cartesian coordinates x—y. 



modeled by two linear springs in the normal and tan- 
gential directions with the spring constants and k^, 
respectively. This harmonic force law approximates the 
interaction for two-dimensional disks for small deforma- 
tions [13, [13 ■ Starting from an arbitrary weakly de- 
formed state, the change in the force exerted by particle 
B on particle A due to an additional small deformation 
of the system is 



= k^ (SI'' ■ n'')rf + fc, [bV ■ 1")^ , 



(3) 



where Sl'^ is the change of the branch vector due to the 
displacement of the particle centers. It is supposed that 
the contacts do not break due to the imposed deforma- 
tion. This is justified by the fact that during the mea- 
surement of the macroscopic elastic moduli, the material 
is subjected only to incremental strain changes so that 
the fabric remains nearly unchanged. The average stress 
increment can be approximated by fl^ [l5| 
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where terms of order (overlap/particle radius) have been 
neglected. This is sufficiently accurate for small pre- 
strains. Here, the sum runs over all contacts Nc within 
the volume V . Using Eqs. ([3]) and @ one obtains 
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A crucial simplification of this expression is obtained, if 
one assumes affine displacements of the particle centers, 
which states that the displacement gradient is constant or 
varies only slowly on the scale of the particle size. In this 
way, the non-affine parts of the relative displacements 
are neglected and the remaining part according to the 
macroscopic strain increment is given by 
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Inserting Eq. (|6]) into Eq. ([5|), one obtains 
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For narrow size distributions and small pre-strain, one 
may replace jZ'^l by the average particle diameter £. Poly- 
dispersity leads, in first order, to an additional factor that 
only depends on the moments of the particle size distri- 
bution J^]. By comparing Eqs. ([7]) and ([2]) the stiffness 
tensor is identified as [13, El 
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Denoting the normal and tangential unit vectors as 
n=(cosa, sina) and sin a, cos a) [see Fig. [TJb)], the 
sum over the contacts located inside the volume V can 
be replaced by an integral over the contact orientation 
distribution P(a) [HI so that 
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where z {=2Nc/N) is the average coordination number, 
(j> {=Ntt£ I is the volume fraction of the packing, and 
P(a) Aa denotes the probability to find a contact with an 
orientation between a and Aa. According to Eq. ([5]), 
the elements of the stiffness tensor are determined by the 
fabric properties: z, and Pip)- 

It is convenient to choose the principal axes of the 
strain tensor increment be as coordinate system so that 
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with (5e^(=Tr(5e=(5y/y) and (^7, the change in the volu- 
metric strain and in the shear deformation, respectively. 
Then Eq. ([2]) can be written as 




/ r C 

c c 

*-'2211 *-'2222 

c c 

^1211 ^1222 

\ c c 

^ ^2111 ^2122 



(5e„-57)/2 



(11) 



The balance of torques requires that the average stress 
tensor remains symmetric, i.e. 5a-^^=5<T^-^^. Hence, the 
following constraints must hold in equilibrium 
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Using © this implies that only those incremental defor- 
mations do not induce torques inside the packing (and 
hence lead to particle rotations), whose principal axes 
are compatible with the contact orientation distribution 
in the sense that 



(sin 2a) = / sin 2a P{a)da = 0. 

J — TT 
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Evaluating © for fabrics that fulfill the constraint 
gives the following expressions for the elements of the 
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stiffness tensor: 
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III. BIAXIAL DEFORMATIONS 

In the following we consider biaxial deformations, 
which means that the principal axes of the strain tensor 
do not change, while the sample is deformed. If the defor- 
mation starts from an isotropic configuration, the fabric 
(P(a)) will remain symmetric with respect to a = 0. 
Hence ([T3|) is fulfilled and (fT7|) vanishes, as one averages 
odd functions with an even distribution. This simplifies 
matters considerably, as it implies that the principal axes 
of stress and strain coincide so that 
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where SP and 5t are the incremental pressure and shear 
stress. 

One defines the bulk modulus E, the shear modulus G 
and the anisotropy modulus A by 
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= E[ ((cos2a)2) + ^((sin2a)2; 



A = ^iiii — ^ = £;(cos2a). 
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For an isotropic fabric, P{a) — the result of Kruyt 
and Rothenburg [3| is reproduced: 

E=An.. G-^(fc„+A:0, A^O. (21) 

ZTT AtT 

In the next section a simple model is proposed for the 
kind of anisotropy, a granular packing develops under 
the influence of shear. 



IV. ANISOTROPY INDUCED BY 
UNILATERALITY 

In order to elucidate the effect of unilaterality on the 
stress-strain relationship, we must consider the closing 
respectively opening of contacts, as a finite volumetric 
strain and shear deformation 7 build up, starting from 
an isotropic, stress-free, jammed packing (unstrained ref- 
erence state with zero overlap). Assuming again that the 
particle displacements may be approximately regarded 
as affine, the distance between the centers of neighboring 
particles changes due to the strain e by 



A^„(a, e)= 
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The change depends on e as well as on the direction of 
the branch vector, a. If the two particles touched each 
other, i.e. ^„ = in the unstrained state, a positive A^„ 
means that the deformation leads to an overlap, while a 
negative value indicates that the particles are no longer 
in contact in the strained state. When there is a gap 
between the particle surfaces in the unstrained configu- 
ration, an overlap can also form, if A^„ is larger than 
the gap. Therefore we extend the notion of an overlap to 
include small negative values, ^„ < 0, which tell the size 
of the gap. 

We introduce the probability density (5(^„,a, e) that 
a particle pair has a branch vector at an angle a with 
respect to the principal axis of the strain tensor e, be- 
longing to the eigenvalue e^, -I- 7, and that the overlap 
respectively the negative gap has a value This prob- 
ability density depends on the strain e. For e = it is 
assumed to be isotropic, i.e. independent of a, and it 
fulfills Qo(Cn) = for ^„ > 0, as there are no overlaps in 
the unstrained configuration, and Qo(Cn) 7^ for ^„ < 0, 
as the unstrained configuration is jammed. We assume 
that the probability distribution simply shifts by A^„, 
Eq. (P^ . under the influence of strain: 



Q(^„,a,e) = Qoi^n - A^„(a,e)) 



(23) 



The probability that a pair of neighbor particles is ac- 
tually in contact is 



/TT pOC 
da / d£,n Q{£,n,a, e) 
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The probability density, that a contact has a certain angle 
a with respect to the principal axis of the strain, which 
belongs to the eigenvalue e„ -I- 7, is 



1 

P'(a,e) = -^y d(.n Q{£.n,a,e). 



(25) 



Applying the assumption (j23p . the integral is in first or- 
der of e given by 

Jo 

= -Qo(0)^(y + |cos2a)(.26) 
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Integrating this over a gives the corresponding first order 
approximation of Af: 



N « -Qo(0)^e.7r. 



(27) 



Hence the properly normahzed first order approximation 
of the probability density of contact directions is 



P(a,e) 



1 
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This approximation can at most be applied for ^ < 1, 

as otherwise the probability density would not be positive 
semi definite, a = is the direction of the principal axis 
belonging to the eigenvalue e^, +7. For e^, < (compres- 
sive strain) and 7 > (by definition) this is the direction 
in which the precompressed system expands during biax- 
ial deformation. Therefore contacts preferentially open 

in this direction so that P{Q, c) ~ 57 ^1 + 7-^ < 57- 

Equation (j28p is the main new result of this paper. In 
this approximation, ((cos2a)^) = ((sin2a)^) = ^ and 
(cos 2a) = 2^ so that the elastic moduli of a granular 
packing are approximately 

Note that due to the presence of the nonzero element 
A in Eq. (fT9| . two independent experimental tests are 
required to determine the elastic moduli of an anisotropic 
material, for example: 



(I) incremental Sey while 6^—0: 



E 



5P_ 
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57=0 



(II) incremental Sj while 5ey=0: 



St 



A^ 



6t_ 

6ev 



5P 
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.57=0 



(31) 



This is in contrast to the isotropic case, where a sin- 
gle experiment with simultaneous incremental Se^ and 
6"f is sufficient to measure both bulk and shear mod- 
uli. In granular media, in contrast to an isotropic elas- 
tic material, a pure shear leads to a pressure increase, 
SP = -ASj > as A < 0. 



V. SIMULATION RESULTS 

We tested the theoretical predictions of SeclIVIbv nu- 
merical simulations. The unstrained initial packing con- 
sists of 3000 rigid disks with particle radii uniformly 
distributed between aniin=0.95 and ainax=l-05 to avoid 
crystalline order. It was generated by a method based 




FIG. 2: Bulk respectively shear modulus (circles respectively 
triangles), determined during an isotropic compression, in 
which the negative volumetric strain — is increased in small 
steps and 7 = 0. Both moduli are given in units of k^/{2n) 
and divided by z(j>- They were calculated from fabric (method 
1, open symbols), respectively from Eq. ((30]) (method 2, full 
symbols). The full line represents z4>, referring to the scale 
on the right. 



on Contact Dynamics simulations, which leads to homo- 
geneous, isotropic, jammed configurations [l^. Such a 
packing is taken as unstrained initial configuration in a 
Molecular Dynamics simulation of soft particles using the 
LAMMPS code [13, [2l|. The interactions between par- 
ticles are modeled with normal and tangential Hookean 
springs (with A;t/fc„=0.5). We checked that the results do 
not depend on the friction coefficient, when it is chosen 
larger than 0.5. The data shown here are for n = 1. 

In a quasi-static compression process, the volume of 
the unstrained packing is gradually decreased by apply- 
ing incremental volumetric strain steps and allowing the 
system to relax between those steps. The process is 
stopped when the total volumetric strain e« — AV/V 
reaches a given value (£„ = —0.04 respectively ey = 
—0.09). This precompressed packing is then sheared in a 
biaxial geometry, while keeping the volume of the system 
constant. The shear deformation is imposed via incre- 
mental steps, and the system is allowed to relax between 
the steps. 

Two different methods are used to determine the elas- 
tic moduli: 

1. Using Eq. ^ the elastic constants can be com- 
puted from the fabric. This formula assumes affine 
deformations and that the contact network remains 
unchanged for incremental strains. 

2. An incremental strain test is simulated by molecu- 
lar dynamics and the elastic moduli are determined 
from Eqs. ([50)1 and ([?T|) . This method allows for a 
change of the contact network and non-affine mo- 
tions of the particles. 

Fig[5] shows that the bulk and shear moduli divided 
by z(f) are approximately constant for strong enough 
compression, — e„ > 0.04, in agreement with Eq. 
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However, as one approaches the unstrained configuration 
(at the jamming transition), the elastic moduh obtained 
from evaluating incremental strain tests soften, whereas 
the ones calculated from the fabric remain unchanged. 
This shows that close to the jamming transition the as- 
sumption that the incremental particle movements are 
afhne fails, in agreement with the findings of [13, [23j . 
Therefore we concentrate on the regime — e„ > 0.04 in 
the following. Remarkably, the ratio of the two moduli, 
G/E is very close to the prediction of the simple theory, 
G/E = (1 + kJkJ/2 = 0.75. Also, the values calculated 
from the fabric agree with the theory, whereas the ones 
obtained by evaluating Eqs.(pO| and ([3T|) are about 20 % 
smaller. 

Figl2] also shows that the contact density {— -i^zcf)) 
increases the more compressed the packing is, while it 
decreases, if a prccompressed sample is sheared at fixed 
volume, see FigUl provided 7 remains small enough. For 
large shear deformation, the system presumably forms a 
shear band and begins to fiow. Then the contact density 
must become independent of 7. The transition from elas- 
tic to plastic response can be seen in FiglH The shear 
stress first increases linearly with 7 and saturates for 
large deformation, indicating the plastic regime. This 
shows, that one can speak about elastic response only 
for ^ < 1. 

e 

Fig|4] confirms that bulk and shear moduli, divided 
by the contact density, do not change, when a prccom- 
pressed state is sheared, in agreement with the theoret- 
ical result Eg.pQ]). The model predicts, however, that 
the anisotropy modulus is proportional to 7/e^ (dashed 
line in FiglS]). Indeed, this is confirmed in the simula- 
tion. The results for A/ E determined from the fabric, 
as well as the ones from the incremental shear tests for 
Eu = —0.04 agree with each other within the error bars 
and show can be fitted by a linear dependence on ^/tv 
. Only the shear test simulation data for e„ = —0.09 
deviate from this behavior for small shear, — 7/et, < 0.7, 
although the evaluation of the fabric perfectly agrees. 
Further simulations are needed to clarify the origin of 
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FIG. 4: Same as Fig[2] but now as function of shear strain 7 
for fixed e„ = -0.09. 




FIG. 5: Ratio of negative anisotropy modulus and bulk 
modulus,— yl/i?, as function of the ratio between shear strain 
and negative volumetric strain, — 7/ei,. Data are shown for 
two volumetric strains, = —0.04 (circles) and — —0.09 
(triangles). They were calculated from fabric (method 1, open 
symbols), respectively from Eg. pO|) (method 2, full symbols). 
The dashed line is the theoretical result, Eq. (|29I) . 



this peculiar behaviour. The simulation data essentially 
confirm the linear dependence of the anisotropy modulus 
on the ratio 7/61,. However, the theory overestimates the 
slope by a factor of about 7. 



T 
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FIG. 3: Shear stress in unites of k^, scaled by as a 

function of shear strain 7 scaled by — e„. Dashed line: — 
-0.04, full line: e„ = -0.09. 



VI. CONCLUSION 

We have shown that the opening and closing of con- 
tacts can explain the anisotropy modulus which is char- 
acteristic for the elastic response of dense granular pack- 
ings under biaxial shear. The theory predicts a linear 
dependence of the anisotropy modulus on the ratio 7/e„, 
which is the only zcroth order combination of the scalar 
invariants 7 — VTre^ — 2 det e and ev — Tre. This is con- 
firmed by simulations, but the theory overestimates the 
anisotropy modulus by a factor of 7. The bulk and shear 
moduli are predicted correctly by the theory, as long as 
one is not too close to the jamming transition, where the 
assumption of affine deformation of the packing fails. 
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